{
 "metadata": {
  "name": "",
  "signature": "sha256:b9bece0d90b43bd9f1b4d084559a7e09e984650f21418289a9991c8b436857b9"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "CDO Tranche Pricing\n",
      "-------------------\n",
      "\n",
      "The par spread $S$ of a CDS tranche is given by\n",
      "$$S = \\frac{\\text{Protection Leg PV}}{\\text{RPV01}},$$\n",
      "\n",
      "where\n",
      "$$\\begin{aligned}\n",
      "\\text{RPV01} &= \\sum_{i} \\Delta_i Z(t_i) \\frac{Q(t_{i-1},K_1,K_2) + Q(t_i, K_1, K_2)}{2} \\\\\n",
      "\\text{Protection Leg PV} &= \\int_{0}^{T} Z(t) (-dQ(t,K_1,K_2))\n",
      "\\end{aligned}$$\n",
      "\n",
      "Here, \n",
      "\n",
      "* $K_1, K_2, T$ is the tranche starting level, final level, and maturity, respectively.\n",
      "* $Z(t)$ is the discount factor curve.\n",
      "* $Q(t,K_1,K_2)$ is the tranche survival curve: $Q(t,K_1,K_2)= 1 - \\frac{f(t, K_2) - f(t, K_1)}{K_2 - K_1}$\n",
      "* $t_i$ is the time of the $i^{th}$ premium payment. And $\\Delta_i$ is the time-fraction of the $i^{th}$\n",
      "  interval, using the Act365 DCF.\n",
      " \n",
      "  \n",
      "And $f(t,K) = E[ \\min (L(t),K)]$ is model-dependent.\n",
      "\n",
      "With the exception for the formula of $Q$, the CDO tranche pricing formula is identical to\n",
      "the CDS pricing formula with $R=0$. So we first define a class for calculating CDS's."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from datetime import timedelta\n",
      "class CDS(object):\n",
      "    \"\"\"A generic CDS class. Can also be used for CDO tranches.\n",
      "    \n",
      "    Attributes:\n",
      "        today (datetime.date): The present date.\n",
      "        maturity (datetime.date): The maturity of the CDO.\n",
      "        remaining_payments (list): A list of datetime.date objects indicating \n",
      "            the remaining payment dates for the premium leg.\n",
      "        R (float): The recovery fraction (0 to 1).\n",
      "        Z (callable function): A function Z(t) that outputs the discount factor\n",
      "            Z(t) for a given time t (datetime.date object). Input and output are\n",
      "            both positive floats.\n",
      "        Q (callable function): A function Q(t) that takes in a single\n",
      "            datetime.date perameter and returns a float representing the tranche\n",
      "            survival probability. This function should be well-defined for all \n",
      "            dates between `today` and `maturity`.\n",
      "    \"\"\"\n",
      "    def __init__(self, today, maturity, remaining_payments, R, Z, Q):\n",
      "        self.today = today\n",
      "        self.maturity = maturity\n",
      "        self.remaining_payments = remaining_payments\n",
      "        self.R = R\n",
      "        self.Z = Z\n",
      "        self.Q = Q\n",
      "    def rpv01(self):\n",
      "        \"\"\"Returns the value of the tranche premium leg for unit notional.\n",
      "        \"\"\"\n",
      "        days = [self.today] + self.remaining_payments\n",
      "        nodes = [(day - self.today).days / 365 for day in days]\n",
      "        qvals = [self.Q(day) for day in days]\n",
      "        total = 0\n",
      "        for i in range(1, len(days)):\n",
      "            delta = nodes[i] - nodes[i - 1]\n",
      "            total += delta * self.Z(days[i]) * (qvals[i] + qvals[i - 1])\n",
      "        return total / 2\n",
      "    def protectionLegPV(self, N = 200):\n",
      "        \"\"\"Returns the value of the protection leg for unit notional.\n",
      "        \n",
      "        Arguments:\n",
      "            N (int, optional): The number of nodes for calculating the integral.\n",
      "        \"\"\"\n",
      "        delta = (self.maturity - self.today).days / N\n",
      "        days = [today + timedelta(days = delta) * n for n in range(N + 1)]\n",
      "        qvals = [self.Q(day) for day in days]\n",
      "        values = [Z(days[i]) * (qvals[i - 1] - qvals[i])\n",
      "                  for i in range(1, len(days))]\n",
      "        return (1 - self.R) * sum(values)\n",
      "    def parSpread(self, N = 200):\n",
      "        \"\"\" Returns the par spread.\n",
      "        \n",
      "        Arguments:\n",
      "            N (int, optional): The number of nodes for calculating the \n",
      "                protection leg integral.\n",
      "        \"\"\"\n",
      "        return self.protectionLegPV(N) / self.rpv01()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Large Homogenous Portfolio (LHP) Model\n",
      "--------------------------------------\n",
      "\n",
      "The LHP tranche pricing model is given by the following equations:\n",
      "$$\\begin{aligned}\n",
      "E[\\min (L(T), K)] &= (1-R) \\Phi_{2}(C(T), -A(K), - \\beta) + K \\Phi(A(K)) \\\\\n",
      "A(K) &= \\frac{1}{\\beta} \\left( C(T) - \\sqrt{1 - \\beta^{2}} \\Phi^{-1} \\left( \\frac{K}{1-R} \\right) \\right) \\\\\n",
      "C(T) &= \\Phi^{-1} (1 - Q(T))\n",
      "\\end{aligned}$$\n",
      "\n",
      "Here:\n",
      "* $\\Phi_{2}$ is the cdf of the bivariate normal distribution. The third perameter being correlation.\n",
      "* $\\beta^{2}$ the correlation\n",
      "* $R$ is the recovery rate.\n",
      "* $Q(T)$ is the survival curve"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import norm, mvn #normal and bivariate normal\n",
      "from numpy       import sqrt\n",
      "\n",
      "def Q_lhp(t, K1, K2, R, beta, Q):\n",
      "    \"\"\"Calculates the Tranche survival curve for the LHP model.\n",
      "    \n",
      "    Args:\n",
      "        T (datetime.date): Should be given in \"days\" (no hours, minutes, etc.)\n",
      "        K1 (float): The starting value of the tranche. Its value should be \n",
      "            between 0 & 1.\n",
      "        K2 (float): The final value of the tranche. \n",
      "        beta (float): Correlation perameter. Its value should be between 0 and 1.\n",
      "        R (float): Recovery rate. Usually between 0 and 1.\n",
      "        Q (callable function): A function Q that takes in a dateime.date and\n",
      "            outputs a float from 0 to 1. It is assumed that Q(0) = 1 and Q is \n",
      "            decreasing. Represents survival curve of each credit.\n",
      "    \"\"\"\n",
      "    if Q(t) == 1:\n",
      "        return 1 # prevents infinity\n",
      "    def emin(K):\n",
      "    # Calculates E[min(L(T), K)] in LHP model\n",
      "        C = norm.ppf(1 - Q(t))\n",
      "        A = (1 / beta) * (C - sqrt(1 - beta * beta) * norm.ppf(K / (1 - R)))\n",
      "        return (1 - R) * mvn.mvndst(upper = [C, -1 * A],\n",
      "                                    lower = [0,0],\n",
      "                                    infin = [0,0], # set lower bounds = -infty\n",
      "                                    correl =  -1 * beta)[1] + K * norm.cdf(A)\n",
      "    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Gaussian Heterogenous Model\n",
      "---------------------------\n",
      "\n",
      "The tranche survival curve is given by \n",
      "\n",
      "$$\\begin{aligned}\n",
      "Q_{\\text{Gauss}} (t, K_1, K_2) & = 1 - \\frac{1}{K_2 - K_1} \\int_{- \\infty}^{+ \\infty} \\phi (z) f(z) dz \\\\\n",
      "f(z) &= \\sigma \\phi \\left(\\frac{\\mu - K_1}{\\sigma} \\right) + (\\mu - K_1) \\Phi \\left( \\frac{\\mu - K_1}{\\sigma} \\right) - (\\mu - K_2) \\Phi \\left( \\frac{\\mu - K_2}{\\sigma} \\right) - \\sigma \\phi \\left( \\frac{\\mu - K2}{\\sigma} \\right) \\\\\n",
      "\\mu &= \\frac{1}{N} \\sum_{i=1}^{N} p_{i}(t,z) F_i (1 - R_i) \\\\\n",
      "\\sigma^{2} &= \\frac{1}{N^2} \\sum_{i=1}^{N} p_i (t,z) (1 - p_i (t,z)) F_i^2 (1-R_i)^{2} \\\\\n",
      "p_{i}(t,z) &= \\Phi \\left( \\frac{C_i (t) - \\beta_i z}{\\sqrt{1 - \\beta_i^2}}\\right) \\\\\n",
      "C_{i} (t) &= \\Phi^{-1} (1 -Q_i (t))\n",
      "\\end{aligned}$$\n",
      "\n",
      "Here:\n",
      "\n",
      "* $K_1, K_2$ are the boundaries of the tranche (they are between 0 and 1).\n",
      "* $N$ is the number of credits\n",
      "* $F_i$ is the fractional face value of the $i$th credit.\n",
      "* $R_i$ is the recovery rate of the $i$th credit.\n",
      "* $\\beta_{i}$ is the $i$th correlation perameter\n",
      "* $Q_i$ is the survival curve of the $i$th credit.\n",
      "\n",
      "Unlike the LHP model, this is heterogenous in that each credit has distinct values of $ R_i, \\beta_i, Q_i$. Moreover, the values of $F_i$ becomes relevant."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import quad\n",
      "def Q_gauss(t, K1, K2, Fs, Rs, betas, Qs):\n",
      "    \"\"\" Calculate the tranche survival probability in the Gaussian heterogenous model.\n",
      "    \n",
      "    Arguments:\n",
      "        t (float): Time. Positive.\n",
      "        K1 (float): Starting tranche value. Between 0 and 1.\n",
      "        K2 (float): Ending tranche value. Between 0 and 1.\n",
      "        Fs (list): List of fractional face values for each credit. Each entry must be\n",
      "            between 0 and 1.\n",
      "        Rs (list): List of recovery rates for each credit. Each entry must be between\n",
      "            0 and 1.\n",
      "        betas (list): Correlation perameters for each credit. Each entry must be between\n",
      "            0 and 1.\n",
      "        Qs (list): Survival curves for each credit. Each entry must be a callable function\n",
      "            that takes in a datetime.date argument and returns a number from 0 to 1.\n",
      "    \"\"\"\n",
      "    Cs = [norm.ppf(1 - q(t)) for q in Qs]\n",
      "    N = len(Fs) \n",
      "    def f(K,z):\n",
      "        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta ** 2)) for C, beta in zip(Cs, betas)]\n",
      "        mu = 1 / N * sum([p * F * (1 - R) for p, F, R in zip(ps, Fs, Rs)])\n",
      "        sigma_squared = 1 / N / N * sum([p * (1 - p) * F ** 2 * (1 - R) ** 2\n",
      "                                         for p, F, R in zip(ps, Fs, Rs)])\n",
      "        sigma = sqrt(sigma_squared)\n",
      "        return -1 * sigma * norm.pdf((mu - K) / sigma) - (mu - K) * norm.cdf((mu - K) / sigma)\n",
      "    emin = lambda K: quad(lambda z: norm.pdf(z) * f(K, z), -10, 10)[0]\n",
      "    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Adjusted Binomial Model\n",
      "-----------------------\n",
      "\n",
      "The adjusted binomial model is given as \n",
      "\n",
      "$$\\begin{aligned}\n",
      "E[min(L(T),K)] &= \\sum_{k=0}^{N} h(k) \\min (Lk, K) \\\\\n",
      "h(k) &= \\int_{- \\infty}^{\\infty} \\phi (z) g(k) dz \\\\\n",
      "g(k) &= \\alpha f(k) \\text{ for } k \\neq l,u \\\\\n",
      "g(l) &= f(l) + (1 - \\alpha)(u - m) \\\\\n",
      "g(u) &= f(u) - (1 - \\alpha)(l - m) \\\\\n",
      "l &= \\text{floor} (m) \\\\\n",
      "m &= p N \\\\\n",
      "u &= l + 1 \\\\\n",
      "p &= \\frac{1}{N} \\sum_{i=1}^N \\frac{(1 - R_i)F_i}{L} p_i \\\\\n",
      "\\alpha &= \\frac{v_E N + (u - m)^2 -((l -m)^2 - (u -m)^2)(u -m)}\n",
      "{v_A N + (u - m)^2 -((l -m)^2 - (u -m)^2)(u -m)} \\\\\n",
      "v_E &= \\frac{1}{N^2} \\sum_{i=1}^{N} \\left( \\frac{(1-R_i)F_i}{L}\\right)^2 p_i(1 - p_i) \\\\\n",
      "v_A &= \\frac{p(1-p)}{N} \\\\\n",
      "f(k) &= \\frac{N!}{(N - k)! k!} p^k (1 - p)^{N - k} \\\\\n",
      "L &= \\frac{1}{N} \\sum_{i=1}^{N} (1 - R_i) F_i\n",
      "\\end{aligned}$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def Q_adjbinom(t, K1, K2, Fs, Rs, betas, Qs):\n",
      "    \"\"\" Calculates the tranche survival probability under the adjusted\n",
      "        binomial model.\n",
      "        \n",
      "        Arguments:\n",
      "            t (datetime.date): Time.\n",
      "            K1 (float): Starting tranche value (0 to 1).\n",
      "            K2 (float): Final tranche value (0 to 1).\n",
      "            Fs (list): List of fractional face values (floats) for each credit.\n",
      "            Rs (list): List of recovery rates (floats) for each credit.\n",
      "            betas (list): List of correlation perameters\n",
      "            Qs (list): List of survival probabilities. These are callable functions that\n",
      "                takes in a single datetime.date argument and returns a float.\n",
      "        Returns:\n",
      "            float: The value of the tranche survival curve.\n",
      "    \"\"\"\n",
      "    if Qs[0](t) == 1:\n",
      "        return 1.0 # initial value -- avoids weird nan return\n",
      "    N = len(Fs)\n",
      "    Cs = [norm.ppf(1 - Q(t)) for Q in Qs]\n",
      "    L = sum([(1 - R) * F for R, F in zip(Rs, Fs)]) / N\n",
      "    def choose(n, k): # Calculates binomial coeffecient: n choose k.\n",
      "        if k == 0 or k == n:\n",
      "            return 1\n",
      "        return choose(n - 1, k - 1) + choose(n - 1, k)\n",
      "    def g(k,z):\n",
      "        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta * beta)) for C, beta in zip(Cs, betas)]\n",
      "        p_avg = sum([(1 - R) * F / L * p for R, F, p in zip(Rs, Fs, ps)]) / N\n",
      "        f = lambda k: choose(N, k) * p_avg ** k * (1 - p_avg) ** (N - k)\n",
      "        vA = p_avg * (1 - p_avg) / N\n",
      "        vE = 1 / N / N * sum([((1 - R) * F / L) ** 2 * p * (1 - p) for R, F, p in zip(Rs, Fs, ps)])\n",
      "        m = p_avg * N\n",
      "        l = int(m)\n",
      "        u = l + 1\n",
      "        o = (u - m) ** 2 + ((l - m) ** 2 - (u - m) ** 2) * (u - m)\n",
      "        alpha = (vE * N + o) / (vA * N + o)\n",
      "        if k == l:\n",
      "            return f(l) + (1 - alpha) * (u - m)\n",
      "        if k == u:\n",
      "            return f(u) - (1 - alpha) * (l - m)\n",
      "        return alpha * f(k)\n",
      "    I = lambda k: quad(lambda z: norm.pdf(z) * g(k, z), -10, 10)[0]\n",
      "    emin = lambda K: sum([I(k) * min(L * k, K) for k in range(0, N + 1)])\n",
      "    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 15
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Exact Value (Recursive, Inhomogenous)\n",
      "-------------------------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def gcd(a, b):\n",
      "    if a * b == 0:\n",
      "        return max(a, b)\n",
      "    if a < b:\n",
      "        return gcd(a, b % a)\n",
      "    return gcd(a % b, b)\n",
      "def Q_exact(t, K1, K2, Fs, Rs, betas, Qs):\n",
      "    Cs = [norm.ppf(1 - Q(t)) for Q in Qs]\n",
      "    N = len(Fs)\n",
      "    g = round(3600 * Fs[0] * (1 - Rs[0]))\n",
      "    for j in range(1, N):\n",
      "        g = gcd(g, round(3600 * Fs[j] * (1 - Rs[j])))\n",
      "    g = g / 3600\n",
      "    ns = [round(F * (1 - R) / g) for F, R in zip(Fs, Rs)]\n",
      "    def f(j, k, z):\n",
      "        if (j, k) == (0, 0):\n",
      "            return 1.0\n",
      "        if j == 0:\n",
      "            return 0.0\n",
      "        ps = [norm.cdf((C - beta * z) / sqrt(1 - beta ** 2)) for C, beta in zip(Cs, betas)]\n",
      "        if k < ns[j - 1]:\n",
      "            return f(j - 1, k, z) * (1 - ps[j - 1])\n",
      "        return f(j - 1, k, z) * (1 - ps[j -1]) + f(j - 1, k - ns[j - 1], z) * ps[j - 1]\n",
      "    I = [quad(lambda z: norm.pdf(z) * f(N, k, z), -12, 12)[0] for k in range(sum(ns))]\n",
      "    emin = lambda K: sum([I[k] * min(k * g, K) for k in range(sum(ns))])\n",
      "    return 1 - (emin(K2) - emin(K1)) / (K2 - K1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 16
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Example (N = 2 credits)\n",
      "---------------------"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from numpy import exp\n",
      "from datetime import date, timedelta\n",
      "N = 2\n",
      "K1 = 0.03\n",
      "K2 = 0.07\n",
      "Fs = [0.5, 0.5]\n",
      "Rs = [0.4, 0.6]\n",
      "def expdecay(today, rate):\n",
      "    return lambda t: exp(-1 * (t - today).days / 365 * rate)\n",
      "today = date(2012,1, 1)\n",
      "Qs = [expdecay(today, 0.0100), expdecay(today, 0.0150)]\n",
      "betas = [0.40, 0.40]\n",
      "tvalues = [today + timedelta(days = 30) * n for n in range(37)] #3 years\n",
      "# LHP perameters average the other perameters\n",
      "R = 0.5\n",
      "Q = expdecay(today, 0.0150)\n",
      "beta = 0.4\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 17
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from time import time\n",
      "x = time()\n",
      "lhpcurve = [Q_lhp(t, K1, K2, R, beta, Q) for t in tvalues]\n",
      "y = time()\n",
      "print(\"lhp: {} ms\".format(round(1000 * (y - x))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "lhp: 119 ms\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = time()\n",
      "gaussiancurve = [Q_gauss(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]\n",
      "y=time()\n",
      "print(\"gauss: {} ms\".format(round(1000 * (y - x))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "gauss: 21154 ms\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x= time()\n",
      "adjustedbinomialcurve = [Q_adjbinom(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]\n",
      "y = time()\n",
      "print(\"adjbinom: {} ms\".format(round(1000 * (y - x))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "adjbinom: 39434 ms\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = time()\n",
      "exactcurve = [Q_exact(t, K1, K2, Fs, Rs, betas, Qs) for t in tvalues]\n",
      "y = time()\n",
      "print(\"exact: {} ms\".format(round(1000 * (y - x))))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "exact: 47635 ms\n"
       ]
      }
     ],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%matplotlib inline\n",
      "import matplotlib.pyplot as plt\n",
      "plt.plot(tvalues, lhpcurve, 'ro', #red\n",
      "         tvalues, adjustedbinomialcurve, 'bo', #blue\n",
      "         tvalues, gaussiancurve, 'go', #green\n",
      "         tvalues, exactcurve, 'ko'\n",
      "         )\n",
      "plt.show"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 22,
       "text": [
        "<function matplotlib.pyplot.show>"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEACAYAAABI5zaHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xt8XNV16PHfMoYQ/JLrxAYJZA2ThAAxtEDAbRwzU3BC\nCrW40Kalo2A1ubek1HbDIzfUjZCUgQCpY5zgS0vygdj54FvShF5MaQMWmPEjNwpP84oxII38kHjc\ngF9yIMZ43T/OmdGZhzRnpHlrfT+f+XjmnH3O7PGR9hytvfbeoqoYY4ypLRPKXQFjjDGFZ427McbU\nIGvcjTGmBlnjbowxNcgad2OMqUHWuBtjTA3K2biLyN0i8qaIPD9Cme+LyKsislVEft+zfZGIvCIi\n20XkykJV2hhjzMj83Ln/CPj8cDtF5AtAUFU/DlwF/Iu7fTpwI/Bp4DygXUSmjbnGxhhjcsrZuKvq\nFmDPCEWagR+7ZX8FTBORWThfCOtVdZ+q7gXWAxeNvcrGGGNyKUTMvQHY5Xm9292Wvr3f3WaMMabI\nCtG4S5bXmmU77nZjjDFFNrEA59gNnOR5fSIw4G4PpW1/PNsJRMQafWOMGQVVzXYj7fvOXch+Jw7w\nIHAlgIjMBfaq6pvAI8ACEZnmdq4ucLeN6OjjjkVVk4/AKbOgLq1QHQROmZUsM+esJlgGdHgey2DO\nWU3JMh2RCIM4fzokHoNARySS8n7ecu3DlLsxFEo5T+JxYzic9/upKn29vXREItwYCtERidDX25tR\nxvtob28fcX+5HpVYr0qsk9Wr+utUKfUaSc47dxH53zh34DNEZKfb5h0DqKr+QFX/S0T+REReAw4C\nf42zc4+IRIGn3LatU52O1eEdLXz3lu+kbHp332TY+2Zqub3w7rGTky8nffx4OKYvtcwxMPnjxydf\ntkajfG3zJg4M7uLNyTBrEKZMPolvRqMph7VGo7R3d9PZ0wPuB2oPBlniKTehoYGDwCTPcQeBCfX1\nyddH+vtT9uOWPzIwkLJtRzzOHQsW0NnTw6TE+3V3s6Sri9mBQEq51W1tHOnvZ+Pbb7Nj0aKU/cYY\n45WzcVfVv/JRZvEw21cDq/1U5OjjjuW7t3yHJUuXpGx/993DWcu/994HyefBmUG6D3U7XzkJh+Dk\nmcHkyyMC6z8p7DwHp9whaHxKWJb298jsQIAlXV0sb2vj8V/8AvnMZ1gSjaY0pN4vgGSDPIovAIDV\nbW3J8+CW7+zpYXlbG+333gtkfgEsA+5YsCDjC8AYY5LK/WeF+6eFDqepaU62CIg2Nc1JlumN9+oJ\n805QTkdpQjkdPWHeCdob702WiSyJKMtQOjyPZWhkSWTY93788ceH3dfX26sdkYjeGA5rRySifb29\nGfuvCwZ1EFRBB0GvCwYzyt0YCqm6ZbyPG8PhZJmOSCR5HgV93D1fRySS8Z4dkYjeGAplrVOxjfT/\nVS6VWCdVq1c+KrFOqpVRL7ftzNquFqJDtajWrLmTCy64hMOH9yW3TZw4jTVr7ky+3rWzn7e6B8Fz\nk//WxEF27ewn0OTc2fbv74cZaSc/Bgb2DzCcUCg07L7ZgUDyznq4/Ym/AI4MDDChvj7jLwAYXYgn\nUStviMdveKeYRvr/KpdKrBNYvfJRiXWCyq1XQsXPLTN//jwee+whmprmUFfXRFPTHB577CHmz5+X\nLLNo0dV8cPhAynEfHD7AokVXJ183TG2AQ2knPwT1U+tJF++L07K0hXBrmJalLcT74qOqe+ILoHPD\nBtrvvTdrI9sajdIeDHLQfZ0I8bRmCfF4pX8BDBfeWd3WNqq6G2Oqm2iOHteSVEJEx1KPurom9u3b\nkXX7nj1Owxzvi7Ng8QJ6zuxJxtyDzwXpWtWVvLvPp1whJTtL3Tv81rQ7/Kx35cFgyl15ezhMZyyW\nce72cJjODRsy36u/nwkNDRnvZYypHiKCDpMKWfFhGT+mT5/Kvn2Z2+vqpiSfB5oC3PSlb/Pf//Zv\neO/Qexx7zLHc9M/fzmiw21a0DTXsAMdAz5k9tK1o497vDx+GGYtChHj8hHcqIXRjjCmR4YLxpXww\nQoeqHxs3btaJE6eldLhOnDhNN27cnFcZVdXQolBqp6v7CC8Kp79tRfHTgZveMavDdMwaY6oDI3So\nVnzM3Q+/cXlvpyzA4cP7UuLy4D82X6i4fKEk7+4jEdrDYZZHIhl35H5z78G5y+9saXHCOi0t7IiX\n9/MZY/JTE2EZcBr4eHzYKefZs2d/1u1796Z2xEavjbL5qs3sPGenJx++kehdQx2cKXH5GU6Z7sXd\nRY3L+5ErvOM3997CN8ZUv5q4c/dj+vSpWbd74/IA6AQOvTAXljfCLXWwvNF5rUP/VSPF5SuZn8wc\nsMwbY2rBuGnc16y5k4kTU9cKSc+XB7jmmhW88fqTcGgn/G4vHNrJG68/yTXXrEiW6d/fnzoaFnLm\nzFcCP6EbyG/qBAvdGFOZaiYsk0siLr9o0dXs3XuAuroprFlzZ0pcHuCXv9wApDdScbq7h9IJk3H5\ntOkOssXl21a00b+/n4apDUSvjZY1bAO5QzdgmTfG1IThelpL+WCM2TKFNHNmU9bpDmbNakqW8TPd\nQW+8V4MXB4emPFiGBi8OppSpVJZ5Y0x1oNazZQpp7twzs24/77yh7cnpDl4C+oCX4K1uZ7qDhGqN\ny0PhM2+MMaU3bsIyfq1ceTtbtz7Lzp07k9saGxtZufL25OuRpjtIZOyMZi6bSlLIzBsbEWtM6Vnj\nniYQCBCLxWhra2NgYID6+nqi0SgBT4PkJ63Sb1y+WvmZ9tji8saUT03MLVNqgcAZ9PW9kLG9qWlO\n8s590+YthL/yBY4cM0jiFnfCock8fvfPmf/Z1E7cSux49SPXnDidLS1cv3Ztxt398kgkZ6euMSa3\nmp9bptT8TEO84rs/5cir03EW13McYTorvvvTlMa9UgdE+ZErdGNxeWPKxzpUR8HPdAdOSuWutCN3\npaRUQnV3vObiZ6riBMuZN6aw7M59lHJNd+C9Yx9pe7V3vI7ET1weLDZvTDHYnXuR+EmpBJg6YRq8\nBdyPs9rs/cBbMGVC9ukSqonfEbE23YExhWd37kXiJ6US4Le76mDjRHjfs0bgyxP57fl1KeWqtdPV\nz4hYi80bU3jWuBeJn5RKgOeffzq1YQd4/zAvvPB08mW8L07oqpAzU6Xb6br5qs3E7opVRQOfi+XM\nG1N4lgpZZrNmBXjrrb4s25t44w2nU7H5y5fy4IfXwWbgADAF+CwsfLeZdfc8UMrqFoWfZQT9lDFm\nvBkpFdIa9zJrbr6UBx9cl7F94cJm1q1zGu4ZZ5zIO7v7YY+nwHSYcWIDv3l+d4lqWlyWM29M/izP\nvYL5ic0f3DEI6YNi98DBD4bLyKk+ljNvTGFZtkyZJWLzkUiEcDhMJBIhFoulxOYnHz0567GTjk5d\naGTT5i0Ezj6Dut9vInD2GWzavKWodS+lfHLmjTEWlqkKfkI3+Ux3UI38xtyt09WMJxZzr3LxeJxQ\nKJQRuvHe4Z/4qVPo738F9noOrIOGhk+w+8XtJa5xceSKy1unqxlvxty4i8hFwEqcMM7dqnpb2v5G\n4B7go8DbQIuqDrj7bgMuBgToUtWvZTm/Ne45xOPxEdMqj5k+mff3pgcu4Oi6SRzaUzux+ZFYp6sZ\nb8bUoSoiE4BVwAXAAPCkiKxT1Zc9xZYDq1X1XhEJAbcCV4rIHwJ/pKqfEhEBfiEi81V10xg/07gT\nCAS4d6T51Q/72x6P76CtbTX9/UdoaJhANNpKIDC7cBUtI+t0NWaInw7Vc4FXVXWHqr4P3Ac0p5U5\nDdgAoKoxz34FjhWRY4EP43yZvFmAeps0n537Rzm3x+M7OP/8m1i7dhux2CbWrt3G+effRDy+o1TV\nLCrrdDVmiJ/GvYHU6Q13u9u8tgKXA4jIZcBkEZmuqt1ADHgd6AceUdXaCABXmB/84C7q61MvS319\nAz/4wV3J19dcs4Jdux4BfoJzWX7Crl2PcM01K0pZ1aJpjUZpDwaTDXwi5t6aNlEZ2CyUpvb5yXPP\nFs9JD5B/HVglIq3AJpyG/LCIBIFPAvXueR4VkUdUtXZy9CpEIBBgy5bNI8blt2x5mGzTEG/5xcMZ\n56vG8E1yojJPp+uSLNkyNgulGQ9ydqiKyFygQ1Uvcl/fgLPi9m3DlJ8EbFPVRhG5HviQqt7s7msD\n3lXV5WnHaHt7e/J1KBQiFAqN/lOZrI6dMpnfDWZ2un5o8iTeOzDU6ZoI3+zadQAnijaLk06awsaN\n36z4Bt4P63g11SoWixGLxZKvOzs7h+1QRVVHfABHAa8Bs3GWlNgKnJpWZgZDXxQ34XwZAHwRWO+e\n42jgUeDiLO+hpvhmNM5QnL+6Uh4zGj+SUq65ealCIK1cQJubl5ap5oV1YyikChmPG8PhclfNmLy4\nbWfWtjtnzF1VPwAWu430S8B9qrpNRDpF5BK3WAjYLiIvAzOBm93tPwN6gReAZ4FnVfU/c72nKY55\noXlQl7axDuaFPpOyyVlFKj0GHc9YRapaWcerGQ9sENM4Eu+Lc/6Xz2fXW7uSo1hPmnkSG+/ZmDJ1\nsJ+ZKqE64/Jgg51M7bARqiYpsejHwP4B6qfWZ130Y8GCz/Hoo10Zx1544QK6utY756nyuHyu0a4p\nZWwqA1OhRmrcc8bcS/HAYu4VpfnKZqUuLTZfhzZf2TxUpsbj8n29vXpdMKiDbjx+EPS6YFD7envL\nXTVjkhhLzN2MP/tkH1wJzAGa3H+vhP0yNO+w37h8PL6DlpZOwuF2Wlo6q2bAlK3raqqdzeduMjRM\nbYDJuMPSXIegXr0djsPNVzNySuWmTTdVRejGpjIw1c7u3E2G6LVRgs8F4ZC74RAEnwsSvXZopOfc\nuWdmPfa884a2OyNiHyN1ROxjVTEi1jJqTLWzDlWTVa6O13g8zrx5n2VgoD+5rb6+gS1bNidHxc6a\nNYe33nox49yzZn2KN954ofgfYgwso8ZUA1tmz+Qt0BTg3u+PMFpT4KhTJsB0kmmVR82ckDZZRe7Q\nTUKlpVXmM5WBZdSYSmR37mZUWpa2sHbKWmfMcsIhiByIJL8U/KwgBdWbVml396bcRrpzt5i7GZX+\n/f2pDTvAMTCwf6jDceXK22lsbEwpkr74N1RvbN4yakwls7CMGZWGqQ1Oh2vanXv91KEOx8Ti3yPN\nVAn5pVVWUujGMmpMJbPG3YxK9Noo3Yu76Tmzx2ngExk1q1LnTs+1gpSjOtMqExk16bNLWkaNqQQW\nljGjEmgK0LWqi8iBCOF4mMiBCF2rujKmMoj3xWlZ2kK4NUzL0hbifZmLYlRrWmU+i4MYU2rWoWqK\nJt4XZ8HiBRl39+lfAvF4nFAoxM6dO5PbGhsbicViFZ9WaXPUmHKyicNMWfjJqEmIx+MjxuardaZK\ny6gxxWR57qYs+vf3O8u4eKVl1CTkis3PnXsmDz7Yl7HdG7qpxLj8cBk1y9vabNUnU1QWczdFk8yo\n8UrLqPHLT1plJcblLaPGlIs17qZo/MxR41cirTISiRAOh4lEIikxechvBalSzVZpc9SYcrGYuykq\nP4uDeMv17++nYWrDsOVGkk9cvlQjYi3mborJOlRNRfObVZOL3+kOLr3071m37j9IvcsP0Nz8pzzw\nwPdG/0GG4SejxpjRsMbdVLR8smpG4ielEqogrdJSJo1Pli1jKlo+WTUj8Tvdgd/ZKkuZVpk1fNPd\nbeEbM2rWuJuy8zNPjV9+pjuoxLRKS5k0hWbZMqbsCplV40clplVayqQpNGvcTdkVcp4aX+9XwLTK\nQqVUWsqkKTTrUDVVoVAZNX75SassZEqlpUya0bDFOkzVa1vRNtSwAxwDPWf20LaiOAtjlHqmyuSy\nfpEI7eEwyyMRa9jNmNidu6kK4dYwsUAsc3s8zIbVmSNQx6pSZ6q0dEnjZamQpuoVMqPGD39plaVN\nqbR0SZMPu3M3VaHUMXc//IyILWRcvrOlhevXrs1Y+Wl5JGLpkuPUmGPuInKRiLwsIq+IyDey7G8U\nkUdF5DkR2SAi9Z59J4nIIyLyaxF5UUQa0483JpdSZ9T4UeqUSkuXNPnIGZYRkQnAKuACYAB4UkTW\nqerLnmLLgdWqeq+IhIBbgSvdfT8Goqq6QUSOA44U8gOY8SPQFBhxOoKUu/sZwCHoXtxdtLt7P6Gb\nfGeqHCl8Y2u2mnzkDMuIyFygXVW/4L6+AVBVvc1T5kXgc6o64L7ep6rTRORU4C5VnZ/jPSwsY8as\nUHPUFFIhZ6q0dEmTbqxhmQZgl+f1bneb11bgcvfNLgMmi8h04BPAPhG5X0SeFpHbRCRrRYwZq/79\n/akNO4xqjppC8pNSCf7CN5YuafLhJ1smW2Ocfpv9dWCViLQCm4B+4LB7/nnA7+N8Qfwb0Ar8KP2E\nHR0dyeehUIhQKOSjasYMKXVGjR8rV97O1q3PZqRUeuPy4D98c4QJvMrH6dcgDUzgSJb7M0uXrF2x\nWIxYLOavsKqO+ADmAg97Xt8AfGOE8pOAne7z84ANnn0twB1ZjlFjxqo33qvBi4PKMpQOlGVo8OKg\n9sZ7y1uv3l6NRCIaDoc1Eolob29mfWbObFKcm6aUx6xZTZ7z9GkweJ3CoIIqDGoweJ329vYly/T1\n9up1waAOOgV0EPS6YFD7srynqX5u25m9LR5uhw41vEcBrwGzce6JtgKnppWZwVD8/iagw30+AXgW\nmOG+vgf42yzvUbL/DFPbeuO9GlkS0fCisEaWRIZt2BPlQotCI5YrlYULm7M27gsXNifLRCIdCi8q\nRBRC7r8vaiTSkSzTEYkkG3b1NPAdkUg5PpYpspEa95xhGVX9QEQWA+vdxvpuVd0mIp3Ak6r6EBAC\nbhGRIzhhmb9zjz0iItcDG9xQ+9PAD3P/PWHM6OTKqIHSZ9X44Sd889prbwPNQI/nyG56ev4k+crS\nJU2CrxGqqvowcEratnbP8/uB+4c59jEge6+SMWUw0jw15cqq8ZNW+eabMVIbdoAe3ngjlnw1oaGB\nl4Bb+Bj9HE8Db/APvGbpkuOQTT9gxp1CrfxUaLkWGjn++En09WXbPjn5/IK/+Spz//VtBo9Mw0mp\nPJt1Ez7Gz//mqwWvr6lsNiukGXeSWTVeZc6q8SMYDA6z/eTk8+Urfsrgke14UyoHj2xn+YqflqKK\npoJY427GnVKv/FQo0Wg0o4EPBoNEo0P1zmdE7JZNWzgj8Dma6hZyRuBzbNm0pQi1NuViE4eZcSne\nF6dtRRsD+weon1pP9Npo1nlq2la00b+/n4apDVnLlFo8Hh8xLu93ROyWTVu45ILvse/wanDHu06b\n2MpDj/098+bPK/rnMIUx0ghVa9yNyaISZ6H0w89MlQBnBD7HC323A7fgjDlsAP6BOU3X8Hx8famq\na8bIGndj8lSJ89T44WeREYCGKRcwMLiD1OybIA1TZrN7/2Olq7AZE1usw5g8VWpGTS7+FhmBg/pr\n4I20o3sYPJK6THehFhoxpWeNuzFZVOI8NX7lSqkEOPkTJ/Dss+mNOwQ/MfT5ss1UuWnTTaNaaMSU\nnmXLGJNFtWbU+HXaaadl3X7qaacmnxdyoRFTehZzN2YY1ZpR40c8HmfBggX09AzF3IPBIF1dXXkv\nAG6hm/KxDlVjiqBaM2oSCpFWWcg1Yk3+rHE3pgiqNaPGr88t+Bxdj3ZlbF9w4QLWdznpkpde+ves\nW/cfpA6cCtDc/Kc88MD3SlPRcWzMC2QbYzJV4spPhXTqpOM4KW3bSe72hHxGxJrSssbdmFGq1jlq\n/Krbt4+NQAQIu/9uBOr27/eUGhzm6Mzt8fgOWlo6CYfbaWnpJB7fUegqGw9LhTRmlKLXRule3J0R\nc4+uqo2MmgkNDcwEvAGmg5AyffDcuWfy4IN9GcemrxFraZVlMNwqHqV8YCsxmSrlZ+WnSlv1yS8/\nS/b19vZqY2NjyupRjY2NGUsJNjcvVQikrTQV0ObmpaX+WDWFEVZisg5VY4qo2jNqkottDwwwob4+\n62LbubJuwNIqi8WyZYwpk1rPqPHL0iqLw7JljCmTWs+oAefuvrOlhfZwmM6WFnbE07NnnNh8Nt7Y\nvI2ILSzrUDWmiKp5jho/dsTj3LFgAZ09Pe6s8NDe3c2Srq6U8I2fBcD9plVa6MYfC8sYU0TVHnPP\npbOlhevXrmWSZ9tBYHkkQnva5GU2IrbwbMpfY8ok0BSga1VX6hw1q7LPP1ON89Qc6e9PadjBWdfp\nyEBm2CnXbJV+0iqHQjdDd/i7dgW45prjbERsGmvcjSmyQFMgZ+dpyh3+DOAQdC/urvg7/AkNDRyE\njDt3by68X4UM3RjrUDWmIrStaBsK3QAcAz1n9tC2oq2s9cqlNRqlPRgkscTHQaA9GKQ1mv9ArsRC\nI5FIhHA4TCQSyVhByu+IWBsNa3fuxlSEal35aXYgwJKuLpZ7cuGXZMmF96sQoRsbDeuwxt2YClDN\nWTWzA4GMztNi8RO6sbi8w7JljKkAtZ5VU0i5s278jYZ1zlXdaZU2QtWYKlDLKz+BZyqD/n4mNDRk\nncqgEPykVEJtpFWO1Lj7ndjrIuBl4BXgG1n2NwKPAs8BG4D6tP1TgN3A94c5f+Fn1DGmxvTGezV4\ncVBZhtKBsgwNXhysionI/ExCVigLFzanTVDmPBYubE4pVwuTmTGWicNEZILbqF8ADABPAn+pqi97\nyvwb8KCq3isiIeDLqnqlZ/9K4CPAO6q6NMt7aK56GDPeVfM8NfkMdhqreDxOKBTKiMunZ97UwmRm\nYx3EdC7wqqrucE92H9CMcyefcBrwNQBVjYnIOs+bnw3MBB4GzhnVJzDGVG1GDeQ32GmsEimVuWaq\n9JNWWc2ZN34a9wZgl+f1bpwG32srcDlwh4hcBkwWkenAXmA50AJcOPbqGjN+VXNGTSEHO/mRK6US\nan9ErJ/GPdstf3oM5evAKhFpBTYB/cBh4GrgP1W1X0SGOxcAHR0dyeehUIhQKOSjasaMH9W88lNr\nNEp7d3fqBGPBIEtGMdipUKpxRGwsFiMWi/kq6yfmPhfoUNWL3Nc34ATxbxum/CRgm6o2isi9wDzg\nCE6n6tHAnaq6LO0Yi7kb40M1Z9T4Wfij1AoxmRnApk1bWLToavbs2c/06VNZs+ZO5s+fV/T6jykV\nUkSOArbjdKi+DjwBXKGq2zxlZuB0lqqI3AQcVtWOtPMsAs62DlVjisfy5QuruflSHnxwXcb2hQub\nWbfuAcBp2C+44BIOH96X3D9x4jQee+yhojfwY1qsQ1U/ABYD64GXgPtUdZuIdIrIJW6xELBdRF7G\n6Ty9uSA1N8bkpVrnqEnws/BHKa1ceTuNjY0p29JDN4sWXZ3SsAMcPryPRYuuLkkdh2ODmIypIeHW\nMLFALHN7PMyG1ZU9c2LWhT+CwYyFP0otV+imrq6JffsyJyarq2tiz57UL6dCp1XafO7GjBPVnFGz\nuq0t2bCDk1nT2dPD8ra2ks1dk02uzJvp06eyb1/m9rq6KSmvS51WaVP+GlNDotdGCT4XdBp4GMqo\nubbyM2pKmQtfSGvW3MnEidNStk2cOI01a+5M2VbqNWKtcTemhiRWfoociBCOh4kciFRNZ2oiF96r\nmLnwhTJ//jwee+whmprmUFfXRFPTnKydqfmsEVuIuegt5m7MOFVpKZOVGnMvlGKsEWsxd2NMikpc\n1q/QC39UmlKPiLU7d2PGoWqehKxa+ZnQLJ+56GGMee7GmNrTv78/tWGHqpmErFoVco3YTZu2EAic\nMeL7WVjGmHGomlMmS7XoRzEUYo3YbCNis7GwjDHjULVOU1Drna5+QjeBwBn09Q2FaEY9t0wpWONu\nTOlV4yRkpVz0o1zyHRFr2TLGmBSBpsCInaeVmFFTrQOd8jHaEbHprEPVGJNVJU5CVq0DnQop24jY\nbKxxN8ZkVYkZNa3RKO3BYLKBT8TcW8u46EepeUfEjsTCMsaYrCoxo6bWBzr5NX/+POLx53FXuMvK\nOlSNMVlVa0bNeDKmlZhKwRp3YypTNWbUjCfWuBtjiqJS7+6reaBTPqxxN8YURSXOUVPrA528bG4Z\nY0xRVGJGzXArOq1uq451ZAvFGndjzKglM2q8ypxRMx4GOvlhjbsxZtQqcVk/G+jksJi7MWZM/GTU\neMsVO6vGYu7uvkpoVK1xN6a2lTqrJpkt4w50smyZMrHG3ZjaVolZNbXAsmWMMWVViVk1tc7mljHG\nFF0lzlMDtT3YycIyxpiiq8SRrLXQ8WphGWNMWQWaAnSt6iJyIEI4HiZyIJK1YY/3xWlZ2kK4NUzL\n0hbiffGi1anWBzv5CsuIyEXASpwvg7tV9ba0/Y3APcBHgbeBFlUdEJEzgX8GpgAfAN9W1X8rYP2N\nMVWi0lZ+qvXBTjnv3EVkArAK+DxwOnCFiHwyrdhyYLWqngl8C7jV3f5b4EuqOgf4ArBSRKYWqvLG\nmNpR6pWfan2wk5+wzLnAq6q6Q1XfB+4DmtPKnAZsAFDVWGK/qr6qqj3u89eBt3Du7o0xJkWpM2pq\nfVUnP2GZBmCX5/VunAbfaytwOXCHiFwGTBaR6aq6J1FARM4Fjk409sYY41XqjJpaX9XJT+OerSc2\nPbXl68AqEWkFNgH9wOHkCUROAH4MfGm4N+no6Eg+D4VChEIhH1UzxtSK6LVRuhd3Z2TURFcV7056\ndiBA+73VM4gqFosRi8V8lc2ZCikic4EOVb3IfX0DoOmdqp7yk4Btqtrovp4CxICbVfXfhznGUiGN\nMbbyU57GNP2AiBwFbAcuAF4HngCuUNVtnjIzgHdUVUXkJuCwqnaIyNHAw8A6Vf3+CO9hjbsxJqdK\nzJcvpzHluavqB8BiYD3wEnCfqm4TkU4RucQtFgK2i8jLwEzgZnf7F4F5QKuIPCsiz4jIGWP7OMaY\n8arUGTXgDHbqbGmhPRyms6WFHfHi5d4Xko1QNcZUjXBrmFgglrk9HmbD6g0Ff79KH8VqI1SNMTWh\n1Cs/VfMoVmvcjTFVo9QrP1XzKFZr3I0xVaPUc9RU8yhWi7kbY2pKITNqqjnmbo27MaamFHrVp0pe\nsm+kxt34G0UCAAAN/0lEQVQW6zDG1JT+/f3OrJJeY5ijptpGsSZYzN0YU1NKnVFTqSwsY4ypKeUY\nxVqu5fos5m6MGVf8zFHjLTeWeWrK2elqjbsxxqQp1B1+Z0sL169dm5IPfxBYHokUPVZvI1SNMSZN\noeapqdSBTta4G2PGpUKt/FSpA52scTfGjEuFyqqp1OX6LOZujBmXCj2StRwDnaxD1Rhjsqj2lZ+s\ncTfGmFGo9JWfLFvGGGNGoRwrPxWKzS1jjDHDKPQ8NaUcyWqNuzHGDCOZUZM2w+Ro5qnJOpK1u7to\nI1ktLGOMMcMo5MpPpV6yzxp3Y4wZRiFXfir1SFYLyxhjzAgCTYERF/lIyaiZARyC7sXdGV8CiZGs\n6XPQFGskq925G2PMGPjNqCn1SFa7czfGmDHwm1EzOxBgSVcXyz0jWZdYtowxxlSmfDJqSrlkn41Q\nNcaYMSjnKFabfsAYY4qolCs/eVnjbowxZeb3Dj+fUaxjbtxF5CJgJU52zd2qelva/kbgHuCjwNtA\ni6oOuPsWAf8IKHCzqv44y/mtcTfG1LSWpS2snbI2IzYfORBJplrmux7rmCYOE5EJwCrg88DpwBUi\n8sm0YsuB1ap6JvAt4Fb32OnAjcCngfOAdhGZlus9jTGm1vhZ+amQo1j95LmfC7yqqjtU9X3gPqA5\nrcxpwAYAVY159n8eWK+q+1R1L7AeuCjvWhpjTJXzs/JTIUex+mncG4Bdnte73W1eW4HLAUTkMmCy\ne9eefmx/lmONMabm+ZmnppDrsfrJc88Wz0kPkH8dWCUircAmnEb8sM9jjTGm5iXmqUnJqlmVmi3T\nGo3ytc2bODC4izcnw6xBmDL5JL45ilGsfhr33UCj5/WJQMrfCKr6OkN37pOAy1X1gIjsBkJpxz6e\n7U06OjqSz0OhEKFQKFsxY4ypWrnmqTkisP6Tws5zSGbUND4lLHNvk2OxGLFYzNd75cyWEZGjgO3A\nBcDrwBPAFaq6zVNmBvCOqqqI3AQcVtUONzTzFHAWTgjoKeBsN/7ufQ/LljHGjHt+Mmq8xpQto6of\nAItxOkNfAu5T1W0i0ikil7jFQsB2EXkZmAnc7B67B4jiNOq/AjrTG3ZjjDEOPxk1fvmaW0ZVHwZO\nSdvW7nl+P3D/MMeuBlbnXTNjjBlnCrnyk035a4wxFcLvyk+JxUFGYtMPGGNMBck1T028L07oqhA7\nz9kJ38bmljHGmFpw6ZebWXfCg07opmP4xt3CMsYYU0Wee7E7s9M1C2vcjTGmikweJHMagyyscTfG\nmCryx5/4QwL3k7OBt5i7McZUkR3xODeFzufA4C5+8o7F3I0xpibMDgT4Zmwjp34hMmI5u3M3xpgq\nNabpB4wxxlQfa9yNMaYGWeNujDE1yBp3Y4ypQda4G2NMDbLG3RhjapA17sYYU4OscTfGmBpkjbsx\nxtQga9yNMaYGWeNujDE1yBp3Y4ypQda4G2NMDbLG3RhjapA17sYYU4OscTfGmBpkjbsxxtQga9yN\nMaYGWeNujDE1yFfjLiIXicjLIvKKiHwjy/6TRGSDiDwjIltF5Avu9okislpEnheRl0TkhkJ/AGOM\nMZlyNu4iMgFYBXweOB24QkQ+mVbsm8BPVPUs4ArgTnf7nwPHqOoZwDnAVSLSWKjKF1ssFit3FbKy\nevlXiXUCq1c+KrFOULn1SvBz534u8Kqq7lDV94H7gOa0MkeAqe7zOqDffa7AJBE5CjgO+B2wf8y1\nLpFKvXhWL/8qsU5g9cpHJdYJKrdeCX4a9wZgl+f1bnebVyfwJRHZBTwELHG3/wz4LfA60AcsV9W9\nY6mwMcaY3Pw07pJlm6a9vgL4kaqeBFwM3OtuPw84DBwPnAxcLyJNo6qpMcYY/1R1xAcwF3jY8/oG\n4BtpZV4EGjyvXwM+ghOrj3i23w38WZb3UHvYwx72sEf+j+Ha7onk9iTwMRGZjRNe+UucO3WvHcCF\nwBoRORU4VlV/IyI7gT8G1orIJJwvitvT30BVs/11YIwxZpTEvXMeuZDIRcD3cMI4d6vqrSLSCTyp\nqg+5DfoPgck4natfV9XH3Ab9R8Bp7qnuUdUVxfggxhhjhvhq3I0xxlSXUY1QFZEDY31jEbnGHdi0\nVUS6ROQkz75F7oCp7SJypWf7TSKyU0T2u6+PiMiatHMdEpH1BajfhSLylIg8JyJPikjYs+8sd2DW\nKyKy0rP9z0TkRRH5wC3zj+7rHhE56JZPOVc56uXZ/tvhzpVHfQp5Ha9y6/+sex03jOb/KK1+Y/3/\nWuU+f84dyLfdrd+zInJpmerkvYYHPfUZVZ0KeQ09+78oIioiG/OtT5b6FfIaPiMinxaRRhE5ICLX\nlqlOZ2U555jqlCFXh+ownaz7R3Nc2jnOx4nNA3wVuM99Ph3oAabh5Mz3ANPcfecCsxLvDxwAnsaJ\n9x8LXATsBPoLUL8zgePd56cDuz37fgWc6z7/L+Dz7vNTgI8DG4BFwC+Aie65PomTNZRyrjLU6yxP\n+cHhzlWm6/h77uuL3PJvlvk6PgVsBSa6208ATnCfHw+8CUwo8zXcn6jDaOtU4Gv4IZzw7PPuz9fG\nCruGv+f+P/0M+Alwbbl/Dz3HjalO6Y9Rzy0jIseJyKOeb6+F7vbZIvJrEfmB+y31sIh8KP14Vd2o\nqu+5L7sZyp3/PLBeVfepkxO/HueXHVV9QlXfTDvVz4Gp7rmuAO7Haehxv6F/ISJPi8gWEfm4u32T\niJzh+SxbRORTafV7TlXfcJ+/BHxIRI4WkeOBKar6hFv0x8ClbrntqvoqTvroR4DfqOph91wvq+ob\n7rkmichG9xv/5yIyy63H4yKy0r0De15EPp3l/22s9fI6AuwXkUeBNcDxIvLf3LqU4zqG3OdXuOfC\nrUu5ruOHgL2qetjd/rqqvu6WPwenwftVma8hwLHuNVyP03j9qVuXclzDi4Eo8BbwRmJHBV3Dd3BS\ntHuAd4CrK+D3EBFpduv0Uvq+0RrLxGHvAZeq6jk4GTHf9ez7GHCHqn4K2AdcnuNcX8H5wYDMQVP9\nZA6aSlCcEbNXuD+0ZwCNOD9YANuAz6rq2UA7cIu7/YfAXwO4P2THqOqLw1VORP4MeFadEboNOAO5\nErIN6gL4JdAozp/y/0tE5rvn+qK7/zJV/TROh/O3Pcd9WFX/APg74J7h6jSGeqV7D+eH71ZgM/Ad\nz75yXMelOBlZFwAvuPvLdR3fwfnCS15DETlXRF4EHgC+UiHXcA5QDzTh/H8s9+wr9TX8Ks7v4AxS\nR6NXyjVcAPxP4CacL6LV5b6GInKcW6dOsn95j4qfVMjhCHCriHwW5w6wXkRmuvviqpr4xXwa54cu\n+0lEWoCzcf40TJw33bC9vqr6ojgDo67ASck8BSfPHpw7qx+7PzTK0Of9GdAmItcDXwZWj1C/03F+\nEBfkWb/3gLOAz+J8+d0nInfg/PAfBXSJiOB8wQ54jvtX93NtFpEpIjJVVTOmbBhDvTJOBfwA55f+\nNcp/Hf8T+CfgA9xfesp3HY8Af4UTZvhjnMbrBuAvcO5w73LfWyjvNXwK50/9C4G73HOX/Bri3HXO\nBTpwGvM/8eyrlGv4f3Aa8kZgJvBVEbmc8v4edgK3q+pvnSahMA38aBt3AVpwvp3/QFWPiEgcNxyC\nM4dMwgee7aknEbkQ+AdgvvutB843XchT7ETg8Rz1eRBYCfwGuAb4H+72KLBBVS8TJ0//cQBVfVdE\nunDuWP8c50/sbPU7Efh34Euq2uep30meYieS+kORpE4gbROwSUT6gRXAVcDfqepnhvks3h8EIcsP\nxljrlXb+JcAlQEhVf1kB1/Gf3OM+Aixzt5f7Oiau4QvAlcAzODHl3wHXq+oz6Yd4357iX8PE7+Jp\nwKM4NzjluIZTcBrJf8IZE/NR4Ig4nYdLqIxr+EWcL5FLcdq/Y4Efquqd6Yd4357iXsPzgMtF5Ds4\n/RwfiMi7WeqUl7GEZaYCb7kNexiY7dmX85tHRP4A+Bdgoaq+7dn1CLBARKaJyHScb8RH0g9P+/dX\nwPs4MULvt+s0hiYx+2tS3Q18H3hCs8x3IyLTcObJuUFVk/FfN8623/3zXHB+2ddl+YizReRjnnNF\ngY3AT4GPishcd99EETnNc9xfuNvn4cQKUzKTClCv9GtzPbDJbdjLfR03AN9yY5h/hNMhB+W7jh8m\n9Rf1fJwO++04naufAvrKfA0F5//nLbeuc9y6ZSubVaGuoXtnewrwNVU9EedO/lful1+lXMP73fc7\nGdgLrFXVO8t5DVV1vqqerKon49ykfnusDXvixPn2Eh8F/D+cjpv/CzyH85/1Es6fOrOB5z3lrwNu\nzHKeLpxv92eAZ4EHPPtagVeBV4ArPdtvw4kBHsb5JXsvy7leBV53t8/F+UV8GvgW0JtWh23AgmE+\n5z/iZAAk6vcM8BF339k48eBXge95jrnUrd+77v/RHpypGV7H+fJ5zj3Xr3Eyaba65/mKe/zjOHf3\niTvDs4tQr9dxYqpHAQfdcw26+34DvFzG67jS/f96xr1mG8p8Hd8DDrnHbwWecN/rGZyf9xcq4Bru\nd+s0iDNJ33rK+LuYdt5ncLNlKuga/oyhrKw7cToxy3YNs5y3nQJly+Q9iElEzgTuUtW5eR1YYUSk\nHqfxSJ+bvmxE5HHgOs38M78Y72XXsQjsGuZvPF/DYsorLCMiVwFrcb61qpaIfAknm2VZrrIllt83\n7SjZdSwqu4Z5GM/XsNhs+gFjjKlBtkC2McbUIGvcjTGmBlnjbowxNcgad2OMqUHWuBtjTA2yxt0Y\nY2rQ/weiurY621YVhQAAAABJRU5ErkJggg==\n",
       "text": [
        "<matplotlib.figure.Figure at 0x7f1259325eb8>"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}